Decoherence of superconducting qubits caused by quasiparticle tunneling 
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In superconducting qubits, the interaction of the qubit degree of freedom with quasiparticles 
defines a fundamental limitation for the qubit coherence. We develop a theory of the pure dephasing 
rate caused by quasiparticles tunneling through a Josephson junction and of the inhomogeneous 
broadening due to changes in the occupations of Andreev states in the junction. To estimate F^, 
we derive a master equation for the qubit dynamics. The tunneling rate of free quasiparticles is 
enhanced by their large density of states at energies close to the superconducting gap. Nevertheless, 
we find that is small compared to the rates determined by extrinsic factors in most of the 
current qubit designs (phase and flux qubits, transmon, fluxonium). The split transmon, in which 
a single junction is replaced by a SQUID loop, represents an exception that could make possible 
the measurement of T^. Fluctuations of the qubit frequency leading to inhomogeneous broadening 
may be caused by the fluctuations in the occupation numbers of the Andreev states associated with 
a phase-biased Josephson junction. This mechanism may be revealed in qubits with small-area 
junctions, since the smallest relative change in frequency it causes is of the order of the inverse 
number of transmission channels in the junction. 

PACS numbers: 74.50. +r, 85.25.Cp 



I. INTRODUCTION 



Over the past several years significant efforts have been 
directed toward designing and implementing supercon- 
ducting circuits with improved coherence properties. For 
quantum computation purposes, the coherence time T 2 
of a qubit must be sufficiently long as to allow for error 
correction^ The unavoidable couplings of the qubit with 
various sources of noise are responsible for decoherence, 
and different types of qubits have different sensitivities 
to a given noise source. For example, the phase and flux 
qubits coherence times are limited by flux noise^^ while 
the transmon parameters are chosen to decrease the ef- 
fect of charge noise in comparison with the Cooper pair 
box.— Flux and charge noise originate from the environ- 
ment surrounding the qubits; in this paper, by contrast, 
we study an intrinsic mechanism of decoherence due to 
the coupling between the qubit and the quasiparticle ex- 
citations in the superconductor the qubit is made of. In 
general one can distinguish two contributions to the time 
T 2 : first, the qubit can lose energy and the correspond- 
ing relaxation time T\ imposes an upper bound to the 
coherence time, T 2 < 27\. Second, additional pure de- 
phasing mechanisms, characterized by the rate T^, can 
shorten T 2 below this upper limit. Recent theoretical^ 
and experimental^— works have highlighted the contri- 
bution of quasiparticle tunneling to the relaxation rate. 
Here we focus on the pure dephasing effect of quasipar- 
ticle tunneling. 

The decoherence rates discussed above are related to 
the power spectral density S(ui) of the noise source: the 
relaxation rate is proportional to the value of the spec- 
tral density at the qubit frequency CJ10, 1/Ti oc S(w\o), 
while the pure dephasing rate is determined by the low- 
frequency part of the spectral density, ~ 5(0) - see, 
e.g., Ref.[Kj. Clearly the latter relationship cannot hold if 



the power spectral density diverges as u — > 0. Because of 
its experimental relevance, a well-studied example of di- 
verging spectral density is that of 1 // noise; in the case of 
1// flux noise, for instance, the decay of the qubit coher- 
ence is not exponential in time, but Gaussian-lik o 10 ' 1 1 (up 
to a logarithmic factor that depends on the measurement 
protocol) . In studying how quasiparticle tunneling affects 
dephasing we find another such example, since the quasi- 
particle current spectral density is logarithmically diver- 
gent at low frequencies when the gaps on the two sides of 
the junction have the same magnitudes (see Sec. IIII|) . We 
show that despite this divergence, a finite dephasing rate 
can be determined. We then estimate the dephasing 
rate for a few different single- and multi-junction qubits 
and find that in most cases is small compared to the 
the quasiparticle induced relaxation rate. An exception 
is the split transmon, in which the two rates can be of 
the same order of magnitude (see Sec. IV A[) . Since it is 
known that quasiparticles limit the relaxation rate in this 
system at sufficiently high temperatures, 9 it may be pos- 
sible to measure the quasiparticle dephasing rate if other 
sources of dephasing can be minimized. 

The quasiparticle dephasing mechanism discussed 
above is due to tunneling of free quasiparticles across the 
junction. Another dephasing mechanism originates from 
quasiparticles weakly bound to a phase-biased junction 
that give rise to subgap Andreev states; the dephasing is 
caused by changes in the occupations of these states that 
make the Josephson coupling and hence qubit frequency 
ui q fluctuate. Because of this additional dephasing, the 
measured decoherence rate 1/T 2 * acquires an inhomoge- 
neous broadening contribution, 1/T 2 * — I/T2, which can 
be suppressed using echo pulse sequences. When the 
average occupation X}L of the Andreev states is small, 
Xq p <C 1, the typical (i.e., root mean square) fluctuation 
of the occupations is given by the square root of x^„. 
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Then for the phase qubit we show in Sec. IVII that the 
typical frequency fluctuation is proportional to the typi- 
cal fluctuation of the occupations divided by the square 
root of the (effective) number of transmission channels N e 

in the junction, ((Aojg) 2 ) 1 / 2 /uj q cx ^Jx^ p /N e . For these 

fluctuations to measurably affect the decoherence rate 
1/T|, the condition {Auj q 2 ) 1 ^ 2 T 2 > 1 should be satisfied; 
using this condition we estimate that this mechanism is 
not a limiting factor to coherence in current experiments 
with phase qubits. On the contrary, it could contribute 
to decoherence in recent transmon experiments,—^ due 
to the small junction area (i.e., smaller N e in comparison 
with phase qubits). However, this possibility will require 
a separate investigation, due to the lack of phase bias in 
the transmon. 

The paper is organized as follows: in the next Section 
we introduce the effective description of a single-junction 
system. In Sec. IIIII we present the master equation 
governing the qubit dynamics and we discuss the self- 
consistent regularization of the logarithmic divergence in 
the dephasing rate. Applications of our results to single- 
and multi-junctions qubits are in Sees. IIVI and [V] re- 
spectively. The role of Andreev states is analyzed in 
Sec. IVII We summarize our work in Sec. IVIII We use 
units h~ks = f throughout the paper. 

II. EFFECTIVE MODEL 

The effective Hamiltonian H for a superconducting 
qubit can be split into two parts, 

H = H + SH, (I) 

where the non-interacting Hamiltonian Hq is the sum of 
qubit and quasiparticle terms, 

Hq = H v + H qp . (2) 

The Hamiltonian for the qubit degree of freedom ac- 
counts for the charging (Ec), Josephson (Ej), and induc- 
tive (El) energies in a system comprising an inductive 
loop shunting a tunnel junction, 

H ip = AE c (N - -Ej cos$+^E L (<p - 27r$ e /$ ) 2 , 

(3) 

with n g the dimensionless gate voltage, <!> e the external 
magnetic flux threading the loop, and $0 = h/2e the flux 
quantum. The operator N — —id/dip counts the num- 
ber of Cooper pairs passed through the junction. The 
quasiparticle Hamiltonian is given by 

H « P = E ^=EE £ ^A' ( 4 ) 

j=L,R 1 = 1 11,0- 

where & na i(6c!<Ti) are annihilation (creation) operators 
for quasiparticles with channel index I and spin a =f, i 



in lead j — L,R to the left or right of the junction. We 
have assumed for simplicity the same number of channels 
N c h and identical densities of states per spin direction uo 
in both leads. Denoting with A 1 the superconducting 

gap, the quasiparticle energies are e 3 n — J (£n) 2 + (A J ) 2 , 

with £J single-particle energy level n in the normal state 
of lead j. The occupation probabilities of these levels are 
given by the distribution functions 

P (H) = ((a$tiKv))<a> = (( &J n\Aii))qp 1 3=L,R, 

(5) 

where double angular brackets ((. . .)) qp denote averaging 
over quasiparticle states. We take the distribution func- 
tions to be independent of spin and equal in the two leads. 
We also assume that 6E, the characteristic energy of the 
quasiparticles above the gap, is small compared to the 
gap, but the distribution function is otherwise generic, 
thus allowing for non-equilibrium conditions. 

The interaction term 8H in Eq. (Q} accounts for tun- 
neling and, as discussed in Appendix A of Ref. 6, is the 
sum of three parts: quasiparticle tunneling Ht, pair tun- 
neling H T , and the Josephson energy counterterm Hej- 
When the superconducting gaps are larger than all other 
energy scales, the only effect of the last two terms is to 
contribute to the renormalization of the qubit frequency^ 
[see also the discussion after Eq. (fl3)) ]; therefore, we ne- 
glect those terms and consider only the quasiparticle tun- 
neling Hamiltonian, 5 H = Ht with 

l,k—l n,m,a 

+ H.c. (6) 

Here the Bogoliubov amplitudes vP n , v 3 n are real quanti- 
ties, since their dependence on the phases of the order pa- 
rameters appears explicitly through the gauge-invariant 
phase difference ip. The elements tik < 1 of the electron 
tunneling matrix t are related to the junction conduc- 
tance by qt — 2gK^2p=i Ip, where gx — e 2 /h is the 
conductance quantum and the transmission probabilities 
Tp (p = I,---, N c h) are the eigenvalues of the matrix 
(2nv ) 2 tP. 

Since we are interested in the dynamics of the qubit 
only, rather than that of a multi-level system, we project 
the Hamiltonian H onto the qubit states |0) and |1), 
which we represent by the vectors (0, l) r and (1, 0) T for 
the ground and excited states, respectively; the two-level 
approximation is justified under the conditions that per- 
mit the operability of the system as a qubit^ (i.e., anhar- 
monicity large compared to linewidth) . Then in terms of 
the Pauli matrices we can write 

K = t^ z - (7) 

where the qubit frequency in general depends on all the 
parameters present in Eq. and, dropping for nota- 



3 



tional simplicity the channel indices^ 

H T =tJ2 [A d nm <r z + A r nm (a + +a-) 

R.c, 



(8) 



+A f I d Lt d fl 

1 nm no ' mo 



where the coefficients A k m7 k = d, r, f , have the struc- 
ture 

A k m = A k c (u^u* - + iA k s (u%u* + . (9) 

Here A k s denote combinations of matrix elements for the 
operators e ±l( ^ 2 associated with the transfer of a single 
charge across the junction, 



Sij — 


(£|sm~|j) 


(10) 


A* = 


2 ( s n _ s oo) 


(11) 


a: = 




(12) 


4 = 


1 / 

2 ( s n + s oo) 


(13) 



and the A k are obtained by replacing sine with cosine 
in the above definitions. As it will become evident in 
the next section, only the terms with k = d and k = r 
contribute to pure dephasing and relaxation of the qubit, 
respectively. 

The term with k = f (in combination with the k = r 
one) contributes to the average frequency shift. More 
precisely, the average frequency shift Slo = Sluej + Sco qp 
has two parts' originating from the quasiparticle renor- 
malization of the Josephson energy and virtual transi- 
tions between qubit states mediated by quasiparticles, 
respectively. The latter part (<5w qp ) is discussed further 
in Appendix [Al Here we note that in the leading (oc i 2 ) 
order, the Josephson part Sluej is the sum of two contri- 
butions with distinct origins. The first one comes from 
the product of the terms proportional to A* nm and A r nm 
in SHt [Eq. ©]■ The second contribution is due to the 
terms we neglected in 5H. (The neglected terms are the 
pair tunneling and Josephson counterterm, as defined in 
Appendix A of Ref. [g.) Since we are studying decoher- 
ence effects in this work, we set A( lm — henceforth. 
Equations Q, ©, and © (with A{ ml = 0) constitute 
the starting point for the derivation of the master equa- 
tion presented in the next section. 



III. QUBIT PHASE RELAXATION: THE 
MASTER EQUATION 

The information on the time evolution of the qubit is 
contained in its density matrix p(t), which we decompose 
as 



p + cj 



(14) 



In this section we present the final form of the master 
equation for the density matrix. The derivation can be 
found in Appendix [Aj where we start from the Hamilto- 
nian of the system presented in the previous section and 
employ the standard Born-Markov and secular (rotating 
wave) approximations^ to arrive at the expressions given 
here. 

The diagonal component p z of the density matrix 
obeys the equation 



dp z 
dt 



= -tr 1 _ f0 + r _vi]p, + [ro-n-ri_M)] (15) 

where, assuming equal gaps in the leads (A L = A R = A), 



2.9T 



e(e 



def(e)(l- f(e + u w )) 
A 2 



V? 



+ 



- AV(e + wio) 2 
e(e + W10) - A 2 



A 2 



\Al 



(16) 



vV-AVte + c^o) 2 



A ; 



r|2 



and ro->i is obtained by the replacement / — > 1 — /. Here 
cjk = e 2 /h is the conductance quantum. The general 
solution to Eq. (fT5|) is 



1^0 



where we introduced the relaxation time T\ as 



-!r- To- 



ll. 



(17) 



(18) 



Equation (|16l) represents the generalization, valid for 
any ujm < 2 A, of the relaxation rate formula derived in 
Refs. li in the limit uj\q <C 2A using Fermi's golden 
rule. Indeed, the assumption that quasiparticles have 
characteristic energies small compared to the gap enables 
us to approximately substitute e — > A in the numerators 
in square brackets in Eq. (|16[) . and neglecting terms of 
order u)\q / A we find 



Tl^O — \A T S \ 2 S'qp(wio) , 



(19) 



where 



S qp (uj) = 
/[(I 



16Ej 



+oo 



dx ■ 



o V%V X + W / A ( 20 ) 



{1 - / [(1 



A) 



and we remind that Ej = Agxf&gK- The agreement 
of Eq. (fill) with the results of Refs. HI validates the 
present approach. Since the relaxation rate is studied 
in detail in those references, we do not consider it here 
any further, except to note that the terms neglected in 
Eq. (|19[) can become important if the matrix element 
is small, \A r s /A r c \ 2 < wio/A. In fact, A r s can vanish at 
particular values of the external parameters used to tune 
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the qubit, for example in the flux qubit when the external 
flux equals half the flux quantum; 5 '- in such a case, one 
needs to retain the term proportional to A r c in Eq. (fT6|) 
to evaluate the (non- vanishing) relaxation rate. 

The master equation for the off-diagonal part of the 
density matrix is 

dp+ „ , 1 



dt 



i (wio + Suj) p 4 



2Ti 



■P+-^4>P+ ( 21 ) 



where Slu is the quasiparticle-induced average frequency 
shifts discussed in the previous section, T\ is defined in 
and the pure dephasing rate is 



Eq. (I 



r. 



4ff T 



A" 
2 



de/(e) [l-/(e)] 



A A 



\J e 2 — (A L ) 2 y / e 2 — (A R ) 2 



+ - 



A A 



(22) 



\/e 2 - (A L ) 2 v / e 2 - (A R ) 2 

where we assumed A R > A L . The general solution to 
Eq. (EH) is 



p+(t)=p + (0)e 



i(vi +6u)t -t/Ta 



(23) 



(24) 



with 

I _ 1 

The pure dephasing rate defined in Eq. ({22} has a 
structure similar to that of the relaxation rate, Eq. (|T6"|) . 
if we substitute wio — > and ^(c) ~~ ^ ^s(c) m tne l at ~ 
ter. Thus we recover the relationship between the power 
spectral density S(ui) of a noise source and the decoher- 
ence rates discussed in the Introduction, Ti^o °c S(u)io) 
and cx S(0). However, in Eq. (|22p we have explic- 
itly assumed an asymmetric junction, A R > A L , and 
extension of this result to the typical case of a symmet- 
ric junction (A R = A L ) is problematic. Indeed, let us 
consider an almost symmetric junction, A R — A L <C A R , 
with \A^\ > \A~\ and a non-degenerate quasiparticle dis- 
tribution [/(e) <C I, e > A R ); then we find, using from 
now on the notation A = A R , 



At 



2\A d A 



dx ■ 



f[(l + x)A] 



^x + {A- A L )/A 



S m (A - A 



(25) 

In the symmetric junction limit A L — > A, diverges 
logarithmically due to the singularity at x = of the in- 
tegrand in Eq. (|25|) ; for example, in thermal equilibrium 
at temperature T ^> A — A L we have 



32Ej 



A1 



-A/T 



In- 



4T 



A - A 1 



IE 



(26) 



Due to the logarithmic divergence, in general we cannot 
simply take oc 5 l qp (0); the correct procedure that leads 
to a finite dephasing rate is presented in the next section. 



A. Self-consistent dephasing rate 

The terms in the right hand sides of the master equa- 
tions (fT5)l and (j2"Tj) are proportional to the square of 
the tunneling amplitude via the tunneling conductance 
gx oc t 2 ; this proportionality is a consequence of the low- 
est order perturbative treatment of the tunneling Hamil- 
tonian [Eq. (jHJ], which enables us to neglect higher or- 
der (in t) terms when evaluating certain correlation func- 
tions involving qubit and quasiparticle operators [see Ap- 
pendix [5] for details]. This implies that those correlation 
functions oscillate but do not decay in time, which is a 
limitation of the used approximation: the inclusion of 
higher order effects introduces decaying factors of the 
from e~ 7 * into the correlation functions, where at lead- 
ing order the decay rate 7 is itself proportional to the 
tunneling conductance. Here we discuss an Ansatz for 7 
whose validity is checked perturbatively in Appendix |B| 
As we show there, a finite decay rate 7 reflects itself into 
a smearing of the singularity for A L = A of the integrand 
inEq. (gSJ), 



+°° dx 

I! X 



+°° dx_ 

'x ./,) 



+00 



dy 



dx 



Vv 

+ °° dy 1 



6(x - y) 



7/A 



^yTT(x-y) 2 + {~f/A) 2 

(27) 

In the problem at hand there are two inverse time 
scales which could serve as a low-energy cut-off to regu- 
larize the integral as in the above equation, the relaxation 
rate Ti-^o and the pure dephasing rate T^. A finite relax- 
ation rate means that the qubit excited level has a finite 
width; one could argue that this uncertainty in the energy 
will in turn reflect itself in an uncertainty of the energy 
exchanged between qubit and quasiparticles, thus smear- 
ing the singularity as in Eq. (|27[) . However, relaxation 
rate and dephasing rate are determined by different ma- 
trix elements [cf. Eqs. (fTTj) - (IT2"1) ]. so one can imagine, at 
least in principle, a limiting situation in which the relax- 
ation rate vanishes, which would then cause the dephas- 
ing rate to diverge. Therefore, we expect that dephasing 
processes will themselves be the ultimate limiting factors 
for coherence, so that 7 = T^. With this identification, 
we arrive at the self-consistent expression for the pure 
dephasing rate 



ME j , Ad , 2 



dx 



+■00 



x{l-/[(l + y)A]}~ 



x)A] 



(28) 



^{x-yY + iT^/Ay 



Equation (|28p is the central result of this paper. It is 
valid for symmetric junctions (or nearly symmetric, A R — 
A L < r ) and we show in Appendix[B]that it agrees with 
the result of the perturbative derivation of the master 
equation extended with logarithmic accuracy to the next 
to leading order in t 2 . 
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Similarly to the relaxation rate, for some specific values 
of the qubit parameters the matrix element A d can be 
small or even vanish exactly. Then one should take into 
account the second term in square brackets in Eq. ([22)) 
to get 

r = J -\A d c \ 2 / dxf[{l + x)A]{l-f[{l + x)A]} 

(29) 

An estimate for the actual dephasing rate is given by 
the larger of the two rates calculated using Eq. (|28l) or 

Eq. d2HD- 



T^/SE, Eq. (H| simplifies to 

r ^K| 2 / +OO ^Re - 1 

7T 1 S| Jo y/x ^Jx + iiy A 

x/[(l + .x)A]{l-/[(l + a; )A]} (33) 
~^K| 2 /o(l-/o)l n ^ 

7T 1 1 1 

We note that both Eqs. ([30| and ([33]) can be written ap- 
proximately in the form 17 oc \A~\ S^T^); however, 
the proportionality coefficients are different in the two 
cases. Solving Eq. (f3"3"|) for by iterations gives 



B. Non-equilibrium quasiparticles 



The relaxation rate in Eq. (fl9j) depends explicitly on 
the qubit properties via the matrix element A r s , while the 
spectral density S qp accounts for the dynamics of quasi- 
particle tunneling. The same structure is present in the 
right hand sides of Eqs. (|28p - ([29|) - a matrix element mul- 
tiplies factors describing the tunneling dynamics. These 
factors can be further simplified under certain assump- 
tions. Here we focus on Eq. (|2"5[) and distinguish two 
cases: first, let us assume that the quasiparticle energy 
is small compared to the dephasing rate, 8E <C T^, and 
that quasiparticles are non-degenerate, /[(l + y)A] <C 1. 
Then integrating first over y and then over x we find 



16E, 



\A1 




(30) 



where 



V2 



dx 



/[(l+x)A] 



(31) 



is the quasiparticle density normalized by the density of 
Cooper pairs. Indicating with /o the typical occupation 
probability, we estimate 1 - x qp ~ foWSE/A. Then solv- 
ing Eq. (f30|) for T^, the requirement T$ » SE can be 
written as 



32Ej 



A d s \ / (l-/o)ln 



ttSE 



8Ej\A d \'f Q (l-f ) 



(34) 

As a specific example, we consider from now on a quasi- 
equilibrium distribution /(e) = e -e / Te , where T e is the 
effective quasiparticle temperature. 18 In this case we have 
SE = T e and /o = e~ A ^ T " < 1, so that the dephasing 
rate is 



W) 



32E, 



,-A/Te 



A 



In- 



7rT„ 



SEj\AfC 



(35) 

To conclude this section, we note that the divergence 
for A R = A L in Eq. (|22|) is a consequence of the square 
root singularity of the BCS density of states at the gap 
edge. Therefore possible modifications of the density of 
states (e.g., broadening 19 ) would in principle lead to dif- 
ferent estimates of the dephasing rate; the effect of a 
small density of subgap states has been recently consid- 
ered in Ref. l20i However, we argue in Appendix [El that 
these potential modifications are not relevant to current 
experiments with Al-based qubits, which we focus on for 
the remainder of the paper. 



IV. PHASE RELAXATION OF 
SINGLE-JUNCTION QUBITS 



-^\A d \ 2 



SE 

/o» x 



(32) 



This condition is in practice difficult to satisfy, since with 
our assumptions fo <C 1, while \A d \ < 1, Ej/A < 1, 
and at the lowest experimental temperatures SE/A ~ 
T/A> 0.01. Thus we conclude that for non-degenerate 
quasiparticles an upper bound for the dephasing rate is 
given by < SE. 

The second case we consider, for both degenerate and 
non-degenerate quasiparticles, is in fact that of small de- 
phasing rate, <§C SE. Then neglecting terms of order 



In this section we consider the dephasing rate for two 
single-junction systems, the phase qubit and the trans- 
mon, under the assumption of small qubit frequency, 
ujiq <§; A (see Appendix IC 1 1 for the flux qubit). The cal- 
culations of the matrix element entering the relaxation 
rate are described in detail in Ref. 0, whose result we 
briefly summarize. Here we use (without giving all the 
details) the same approach of that work to obtain the 
matrix elements for dephasing. Interestingly, in all cases 
the pure dephasing rate T$ turns out to add at most a 
small correction to I/T2 in comparison with the relax- 
ation term I / 2Ti . 
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A. Phase qubit 



B. Transmon 



In a phase qubit, the charging energy Ec is small com- 
pared to the transition frequency wio- The latter depends 
on the external flux via the position ipo of a minimum in 
the potential energy of the Hamiltonian in Eq. ([3]), as 
determined by 

Ej sin ipo + E L (ip - 2^$ e /$ ) = (36) 

Then the frequency is 

^8E C {E L + Ej coscpo). (37) 

For a small effective temperature T e <C the relaxation 
time is 



f 



1 w, 



Ti 7T Wio 



P_ e -A/T e 



7rT e 



(f + cost^o) 



where 



VSE c Ej 



(38) 



(39) 



is the plasma frequency of the junction. 

Within the same approximations used to obtain the 
above formulas^ the matrix element for dephasing is 



Ec 



(1 - costpo) 



(40) 



and substituting into Eq. (|35|) we get 



E c w 
2ir to 



P_ p -&/T c 



A 

tZ 



In- 



87rT e w 2 



E c up{l - cosv?o) 



(41) 



x(l — COSf ). 



Note that the factor in front of e~ A / Tc is smaller for 
in comparison with that for 1/Ti because the matrix el- 
ement for dephasing is smaller than that for relaxation 
by a factor Ec /wio- At low temperatures the terms in 
square brackets in Eq. (|4"Tj) are dominated by A/T e and 
hence, neglecting factors cos ip as they are small com- 
pared to unity, the condition 2TiT,p > 1 can be written 
as 



Te < fwioy/ 3 (Ec_ 



V A J 



^10 



2/3 



(42) 



Typically for a phase qubit the product on the right is 
of order 10 -2 , while T e /A ~ 10 — 1 . Therefore the pure 
dephasing contribution to Ti [Eq. (|24|) ] can be neglected. 
Interestingly, for a quasiparticle temperature of the order 
of the base temperature, T/A ~ 10~ 2 , relaxation and 
pure dephasing would have similar order of magnitudes, 
although both would be much smaller than at T e /A ~ 
10 _1 due to their common exponential suppression by 
the Boltzmann factor. 



The Hamiltonian of transmon is given by Eq. © with 
El = 0, supplemented by a periodic boundary condition 
in phase. 4 For our purposes, the transmon can be consid- 
ered as a particular case of the phase qubit with ipo = 
[see Eq. (|36p]. With these parameters, one obtains from 
Eq. (|38|) the correct estimate for the relaxation time Xi , 



1 

T7 



-0J v e 



-A/T c 




(43) 



However, the vanishing for ip Q — of the matrix clement 
in Eq. ([4T)]) is not the correct result for the transmon: 
careful evaluation of the matrix element, following the 
procedure outlined in Appendices B and C of Ref. [f| gives 

an exponentially small value, oc exp — yj 8Ej / Ec ■ 

This exponential suppression is sufficient to ensure that 
the dephasing rate is dominated by the contribution in 
Eq. (|29[) , since the matrix element entering that equation 
has no such suppression, 



l4f| a =T 



Ec 



l_Ec 
32 Ej 



(44) 



Substituting this expression into Eq. (|29|) . for the quasi- 
equilibrium distribution function we find 



7T ZA 



(45) 



Using Eqs. (|43j) and (|45l) it is easy to show that for the 
transmon 2TiT0 1; therefore, as for the phase qubit, 
the pure dephasing contribution to T2 is negligible. 



V. PHASE RELAXATION OF 
MULTI-JUNCTION QUBITS 

The results of Sec. HH are readily generalized to multi- 
junction systems by following the same procedure as in 
Sec. V of Ref. @. Assuming the same gaps and distribu- 
tion functions in all superconducting elements, we simply 
need to substitute 



Ej 



A 



M 

E 

3=0 



E 



J3 



A 



(46) 



in Eqs. (|28l) and (|29|) (and hence in subsequent equations 
in Sec. IIII Bp . Here index j denotes the M + 1 junctions 
with Josephson energy Ejj and capacitance Cj , while the 
matrix elements are defined by 



<^i((l|sin||l)- 



|sin^|0) j ,47. 



with tfj the phase difference across junction j. The sim- 
ilar definition for Aj? • is obtained by replacing sine with 
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cosine. We remind that the phases are not independent, 
as they are constrained by the flux quantization condition 



M 

j=o 



27T/, / = $e/$0 



(48) 



Below we consider explicitly the two-junction split trans- 
mon, while the many-junction fluxonium is analyzed in 
Appendix IC2l 



A. Split transmon 

The split transmon single degree of freedom is governed 
by the same Hamiltonian of the single-junction trans- 
mon, but the SQUID loop has a flux-dependent effective 
Josephson energy 

Ej(f) = (Ejo + Ej X ) cos (tt/) ^l + ^tan 2 (irf) (49) 



with 



d = 



\E 



JO 



E 



Ejo + E, 



71 



(50) 



quantifying the junction asymmetry. In quasi- 
equilibrium at the effective temperature T e , the relax- 
ation time is given by£ 



Ti(f) 



-A/T, "K/)+"p(°) 



where 



cu p (f) = y/8E c Ej(f), E C = 



2(C + d) 



(51) 



(52) 



We note that the smaller the asymmetry, the larger the 
tunability of the qubit, since w p (0)/w p (l/2) = \j\fd. 
However, this flexibility comes at the price of enhancing 
the relaxation rate, T x (0)/Tx (1/2) = (1 + d)/(2d 3 / 4 ). In 
Fig.[T]we plot the normalized relaxation rate Ti(0) 
as a function of reduced flux / for three values of the 
asymmetry parameter. We note that the relaxation rate 
rises by about a factor 1.5 up to / ~ 0.4, but can increase 
sharply for small asymmetry as / — > 0.5. 

The matrix elements for dephasing are [cf. Eq. (|40[)] 



\A d \ 2 -- 



Eg 



[1 - cos(tt/ ± ■&)} 



(53) 



where the upper (lower) sign should be used for j = 1 
(j = 0) and 



tan($) = dtan(7r/) 



(54) 



Note that in contrast with the single junction transmon, 
the matrix elements in general do not vanish (except at 
/ = 0). Using Eqs. (g9j, ([53]), and (JM) we find 



E- 

i=o 



Ej 3 \A a a 



64 



:Er 



<4(o) 



1 



(55) 



fC 4 




/ 



FIG. 1: Normalized relaxation rate Ti(0)/Ti(/) vs. reduced 
flux / for (top to bottom) d = 0.02, 0.05, 0.1. Inset: nor- 
malized frequency uj p (f) /uj p (0) vs. reduced flux for the same 
values of the asymmetry parameter (but decreasing top to 
bottom). 



.- l.o - 




FIG. 2: Normalized dephasing rate 2Tir^ vs. reduced flux 
/ for (top to bottom) d — 0.02, 0.05, 0.1. Other parameters 
are specified in the text after Eq. (|56|) . The vanishing of 
as / — > is an artifact of the approximations used to obtain 
Eq. (|56|) ; a finite dephasing rate at any flux would be obtained 
by including a subleading contribution analogous to Eq. (|45|) . 



and the above-described generalization to multi-junction 
systems of Eq. (|35|) gives 



r ^ = 2^ EC 



,-A/T e 



+ ln 



E c (^(0)/^(f) - 1) 



(56) 



In Fig. [2] we show examples of the dependence of 2Tir^ 
on flux for different values of the asymmetry parameter 
d and typical values of the other dimensionless parame- 
ters (Ej(0)/E c = 80, w p (0)/A = 0.2, T e /A = 0.06); we 
note that near / = 1/2 and for small asymmetry, pure 
dephasing dominates over relaxation, 2TiT ^ > 1. There- 
fore the pure dephasing effect of quasiparticle tunneling 
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could be measured in a split transmon if other sources of 
dephasing (such as flux, photon, and charge noise) can 
be suppressed. Charge noise, in particular, can become 
the dominant dephasing mechanism as / — > 1/2, since the 
Cooper pair box regime of small Ej(f) / Ec is approached 
in this case for small asymmetry. 4 However, the contribu- 
tion of to 1 /T% becomes relevant and thus potentially 
observable at values of reduced flux smaller than 1/2, 
where the system is still in the transmon regime; for ex- 
ample, for / - 0.35 where Ej(f)/E c ~ 0A5Ej(0)/E c , 
we estimate 2T 1 Ta ~ 0.4. 



VI. T 2 * AND ANDREEV STATES IN A 
JOSEPHSON JUNCTION 



ture, in Eq. (j37j) we replace Ej by 



A Nch 

P =i 



1 - 2xi 



(58) 



From this substitution we see that a change in the oc- 
cupation of a single Andreev level can lead to a small 
change 8Ej in the Josephson energy and hence in the 
qubit frenquency, with a relative frequency shift of the 
order of 5Ej/Ej ~ 1/N c h. This effect could be mea- 
surable in small junction (N c h ~ 10 5 ) qubits and may 
have already been observed in a transmon, where slow 
frequency jumps of few parts per million magnitude have 
been measured^ More generally we find for the qubit 
frequency uj q at a given set of occupation numbers x£ 



In the previous sections we have considered the pure 
dephasing due to the interaction between tunneling 
quasiparticles and qubit. Here we study a different 
quasiparticle mechanism affecting the measured dephas- 
ing rate 1/T 2 * : as discussed briefly in Sec. Inland in more 
detail in Ref. [|| the quasiparticles renormalize the qubit 
frequency by shifting it by an amount 5u> which depends 
on the quasiparticle occupation. Therefore fluctuations 
in the occupation induce frequency fluctuations that can 
cause additional dephasing. In this section we focus on 
the phase qubit and show that this mechanism is not 
active during a single measurement, so that it does not 
contribute to the pure dephasing rate r^; however, it 
can contribute to the time T 2 * by changing the qubit 
frequency from measurement to measurement. In other 
words, this mechanism being slow on the scale of the 
qubit coherence time, its dephasing effect can be cor- 
rected by using echo techniques. 

In a Josephson junction, weakly bound quasiparticles 
occupy the Andreev states that carry the dissipation- 
less supercurrenti22 Changes in the occupations of these 
states affect the value of the critical current (or equiva- 
lently of the Josephson energy) and in turn fluctuations 
in Ej lead to frequency fluctuations. As we show be- 
low, the parameter determining the relative magnitude 
of these fluctuations is the inverse square root of the 
(effective) number of transmission channels through the 
junction; therefore this fluctuation mechanism could be 
relevant in small junctions. For each transmission chan- 
nel p (p = 1, ■ • ■ , N c h) with transmission probability T p 
[defined after Eq. ([5])]. we find a corresponding Andreev 
bound state with binding energy [see Appendix ID] 



= A - Ei 



E„ = A 1 - 



1 

-T p sin" 



(57) 



This result is valid for T p < 1; the expression valid for 
arbitrary T p can be found in Ref. l22l The (zero tem- 
perature) Josephson energy entering Eq. © is given by 
Ej = A J2 P Tp/4. To account for the occupations x£ of 
the Andreev states, due for example to finite tempera- 



wio 



&Ec A A 

COSlfio > — L V X„ 



^ 4 

p=l 



(59) 



Here we assumed that on average the occupation num- 
bers are small, x^ p = (Xp) <C 1; in quasi-equilibrium 
the average takes the exponentially small value x^ p = 
e -A/T Sj 23 p rom this expression we see that fluctuations 
of the occupations of the Andreev states lead to frequency 
fluctuations. The mean square fluctuations of x£ are re- 
lated to the average x^ as^ 



((Ax^Y) ee ((xf 



(60) 



Using this expression for the non-degenerate case x^ p <C 
1, we find for the root-mean-square frequency fluctua- 
tions 




(61) 



where 



T T 2 



(62) 



is the effective number of channels; N e coincides with 
N c h if all the channels have equal transmission probabili- 
ties. The number N e can be estimated independently by 
measuring the so called subgap structure due to Andreev 
reflections^ 



SIi gr 
Sh 2g K 



(63) 



where the first factor in the right hand side is the ra- 
tio between the current step 5I\ measured as the volt- 
age increases from below to above 2A/e and the sub- 
gap current step SI 2 at V ~ A/e. This ratio is re- 
lated to junction transparency and is of the orde r 26 i 27 
SI2/SI1 ~ 10 -5 — 10~ 3 , while depending on junction area 
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the ratio between junction conductance gx and the con- 
ductance quantum gx is grl 9k ~ 1 — 100, so we estimate 
iV e ~ 10 3 to 10 7 for junction sizes from small to large. 

The dephasing effect of the above frequency fluctua- 
tions gives observable contribution to T 2 * if 



((A0 2 > 1/2 ?2 > 1 
Using Eq. (IBTj) this condition is 



(64) 




oit~nb- ( 65 ) 



Assuming equilibrium between the occupation factors of 
Andreev states and free-quasiparticle states at the effec- 



tive temperature T e s» 140 mK (so that x, 



A _ P -A/T c 

qp — 



since w m 



10 11 s" 1 we find T 2 > 10" 6 s (10" 4 s) for 



small (large) juctions. For phase qubits, which are fabri- 
cated with large junctions, this estimate is two to three 
orders of magnitude longer than the observed coherence 
time^ Therefore fluctuations in the occupations of An- 
dreev levels do not contribute significantly to dephasing 
in current experiments with phase qubits. 

The dephasing effect of the frequency fluctuations can 
be corrected using an echo pulse if the occupations do not 
change during a single measurement. In other words, if 
the rate at which the occupations change is small com- 
pared to I/T2, then the fluctuations contribute to the de- 
coherence time T 2 * rather than to T 2 . Within our model 
Hamiltonian, Eq. (p}, the only processes that can change 
the quasiparticle occupations are due to the interaction 
between qubit and quasiparticles; for an occupied An- 
dreev level, this interaction leads to its ionization, with 
the qubit relaxing and giving its energy to a bound quasi- 
particle which is then excited into the continuum part of 
the spectrum. Since this process relaxes the qubit, it 
can in principle contribute to 1/T\. We show in Ap- 
pendix [D] that this intrinsic contribution is small com- 
pared to the relaxation rate due to the interaction of the 
qubit with the bulk quasiparticles. There are of course 
extrinsic mechanisms that could affect the occupations 
of the Andreev states and hence the rate of frequency 
fluctuations. An example of such a mechanism is flux 
noise; we estimate that the ionization rate due to flux 
noise is in fact small compared to the experimental 1 /T 2 
- see Appendix ID 21 Another mechanism is the quasipar- 
ticle recombination caused by the electron-phonon inter- 
action. The recombination rate is ~ x qp /rQ, with the 
characteristic time To ~ 10~ 7 — 10~ 6 s in aluminum and 
~ 10 -10 s in niobium.] 28 ' 29 Since at low temperatures 7 -^ 
x qp ~ 10~ 7 — 10~ 8 , we find that the recombination rate 
is much smaller than 1/T 2 . 

So far we have considered the effect of fluctuations 
of the Andreev levels occupations. Other mechanisms 
can in principle contribute to decoherence. For exam- 
ple, fluctuations of the order parameter A in the vicin- 
ity of the junction also affect the Josephson energy, see 
Eq. (|58p ; however, at low temperatures the typical time 



scale over which A changes in response to a sudden per- 
turbation is very short, of order 1/A|2! so these fluctu- 
ations do not lead to additional decoherence. Another 
mechanism is associated with fluctuations in the number 
of free (rather than bound) quasiparticles. As discussed 
at the end of Sec. [TIJ there are two contributions to the 
average frequency shift - the Josephson one, 5u>Ej, and 
the quasiparticle one, Suj qp . Fluctuations of free quasi- 
particle occupations affect the latter, but their contribu- 
tion to inhomogeneous broadening is small. Indeed, the 
average frequency shift can be obtained by considering 
the effect of quasiparticles on the junction impedance^^ 
in quasiequilibrium the contribution of the normalized 
quasiparticle density a; qp to the quasiparticle part Y qp of 
the junction impedance Yj is smaller than the term in Yj 
proportional to x qp by the parameter \jT e jixi\Q. More- 
over, the root mean square fluctuations of a; qp scale as the 
inverse square root of the volume of the electrodes^ 4 - and 
can therefore be neglected for macroscopic electrodes. 



VII. SUMMARY 

In this work we have studied decoherence caused by 
quasiparticles in superconducting qubits and obtained 
estimates for the pure dephasing rate Ta, and for the 
contribution of inhomogeneous broadening to the deco- 
herence rate 1/T 2 . We have presented a master equa- 
tion approach that not only reproduces and generalizes 
the formula for the relaxation rate 1/Xi of Refs. 00 [see 
Eq. ([175)) )]. but also gives a self-consistent expression for 
the pure dephasing rate T^, Eq. (|28j) . Moreover, in study- 
ing 1/T 2 we have derived a formula, Eq. dTJTI) . for the typ- 
ical fluctuation of the qubit frequency due to change in 
the occupations of Andreev states. These two equations 
are our main results. 



Application of Eq. (|2"5)l to single-junction qubits such 
as the phase qubit, the transmon fSec. IIVI) . and the flux 
qubit (Appendix IC 1[) . and to the many-junctions flux- 
onium (Appendix IC 2|) shows that in these systems the 
pure dephasing rate is a small contribution to decoher- 
ence, 2TiF0 < 1. In the split transmon (Sec. IV Al) . on 
the other hand, the quasiparticle dephasing rate can be 
larger than the relaxation rate when the external flux 
that tunes the qubit frequency approaches half the flux 
quantum, see Fig. [5] together with its temperature and 
flux dependence [Eq. (|56|) ]. the increased importance of 
in this regime could permit its experimental measure- 
ment. 

Finally in Sec. I VII we have considered the contribu- 
tion to the decoherence rate 1/T 2 due to quasiparticles 
bound into Andreev states localized near the Josephson 
junction. Fluctuations of the occupations of these lev- 
els from measurement to measurement can in principle 
induce dephasing which can be corrected with an echo 
pulse. In practice, this mechanism gives negligible contri- 
butions to dephasing in current experiments with phase 
qubits: due to the short observed T 2 time, Eq. (f6"4")) im- 
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plies that the fluctuations of the occupations would need 
to cause relative frequency fluctuations of the order 10 -3 
to start affecting the coherence of the qubit. 
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time evolution, such as 



(iwio + e^-dT 



+ t{±At nP ±(t) [tf (1 - /*) + (1 - ft;)!*] \AZ n 

x [(1 ± Pz (t)) /«) -(1 T p,{t))0- am 

(A5) 



The terms in curly brackets originate from averages of 
one qubit operator times four quasiparticle operators 
evaluated in the Born approximation^ for example 



Appendix A: Derivation of the master equation 

In this Appendix we summarize the main steps of the 
derivation of the master equation presented in Sec. Mil 
Our starting point is the von Neumann equation*" which 
we write for the two components of the qubit (i.e., re- 
duced) density matrix as 



. ry ■ n ' n 



hj,p 



p+(t) -/£) + (! 



(A6) 



where {•; •} is the anticommutator. The solution of 
Eq. (|A5]l is 



^ = -iTr{[6H;p t ]a*} (Al) 
- iuj wP+ - iTr { [5H; p t ] o+ } (A2) 



Here pt is the total density matrix of the system, com- 
prising both qubit and quasiparticles, [•; •] denotes the 
commutator and, as discussed in Sec. El for our pur- 
poses the interaction Hamiltonian SH — Ht is given by 
Eq. (JSJ) with A{ lm = 0. More useful forms of the traces 
in the right hand sides of the above equations are 



Tr- 



Ht\ Pt 



2t«(<r 



and similarly 



a L ^a R ) 
tn n<j ma I 



H.c.' 



(A3) 



(A4) 



where angular brackets denote quantum statistical aver- 
aging with respect to the total density matrix and the 
prime denotes that Hermitian conjugation is not applied 
to qubit operators (i.e., Pauli matrices). 

The averages in the right hand sides of Eqs. (|A3|I - (|A4[) 
can be found by solving the equations governing their 



it dre K±"io+^-^+iO+)(t-r) 

Jo 



± KnP±(j) [/„ L (1 - C) + (1 - ti)f*" 



1 



A r * 



[(1 ± p,(r)) / n L (l - /*) - (IT Pz (r))(l 1 



(A7) 



A similar expression can be derived for the average in 
Eq. (j A4|) that contains a z . After substituting these 
expressions into Eqs. (|A3j) - (|A4j) and the results into 
Eqs. (|A1|) - (|A2|) . we perform two additional approxima- 
tions. First, we neglect fast rotating terms; this so-called 
secular (or rotating wave) approximation 15 is valid when 
the decoherence rate is small on the scale of the qubit fre- 
quency, I/T2W10 <C 1, and it amounts to keeping in the 
equation for p z only the terms proportional to (1 ± p z ) 
and in the equation for p + only those proportional to p+ . 
With this approximation we find 



d Pz (t) 
dt 



-2? fdr^K 

JO r> m 



e + e 



(A8) 



+ [/n(l-/^)-(l-/ I f)/™ 1 



11 



and 

dp+{t) 
dt 



ius w p+(t) 



(A9) 



2 



x{2|^ m | 2 [e++- + 



3 +— 



where we use the shorthand notation 



(A10) 



Next we introduce the Markov approximation 15 by 
substituting in the integrands of Eqs. (| A8|) - (| A9|) p z (r) —> 
p z {t), P+{t) — > e - tu >io(t-T) p + (t} anc j extending the lower 
integration limits from to —00. Then the r-integrals 
can be performed using the identity 



dTe i(u,+iO+)(t-T) =i p± +7T S(u;) 



(All) 



where P denotes the principal part. We note that in 
Eq. (|A8[) the contributions of the principal parts cancel 
out, while after rewriting the summations over n, mas 
integrals over the quasiparticle energies the <5-functions 
can be used to eliminate one of these integrals. Assuming 
equal gaps in the leads, we finally arrive at Eq. (|15[) . 

Applying the above steps to Eq. (IA9|) . we find that the 
principal parts cancel out in the term proportional to 
A^ m ; in that term we assume different gaps with A# > 
Al to get expression (f2"2")l for the pure dephasing rate T^. 
On the other hand, we can take the gaps to be the same 
in the term proportional to A r nm ; then the (5-functions 
give rise to the contribution — l/27ip + in Eq. pTj) . As 
for the principal parts, they contribute a term idujp+(t) 
with 



5Q = \A r f [F qp (-cj 10 ) - F qp (oj w )} 



(A12) 



The function F qp is defined in Appendix A of Ref. |6|; as 
in that work, we have neglected here contributions sup- 
pressed by the factor wio/A. We note that while 8ui has a 
structure similar to that of <5w qp in Ref. [H, due to the pro- 
jection onto the qubit subspace described in Sec. [TT] the 
expression in Eq. (IA12[) accounts for virtual transitions 
between the qubit states only and neglects those to other 
states of the full system. In systems with small anhar- 
monicity (e.g., the transmon and phase qubit) these tran- 
sitions cannot be neglected and the average frequency 
shift must be calculated using the formulas in Ref. [g. Fi- 
nally, we remind that the total average frequency shift 
8to contains also a Josephson part Sue,, as discussed in 
Sec. Ill 



t) perturbative considerations of Appendix [Al in order to 
regularize the logarithmic divergence in Eq. (|22[) for equal 
gaps. Here we focus on the next to leading order contri- 
butions to validate that equation. First, however, let us 
discuss briefly the smearing of the singularity, Eq. (|27[) . 
which is obtained as follows: after the Markov approx- 
imation, the term in Eq. (|A9[) proportional to Af lm is 
explicitly 



iP P+ (t) [fn (1 ~ O + (1 " /»)/*] Km\ 
r e i(4-£m+n)(t-r) + e i(-e£+e£-M7)(t-T) 



lim / dr 

7^0+ 



(Bl) 



Rather than taking the limit, we assume 7 small but 
finite (in particular, 7 <C ujiq for the rotating wave ap- 
proximation to be valid). After integration the last line 
becomes 



27 



(e L 



)2 +7 2 



(B2) 



This explains the origin of the last factor in the second 
line of Eq. (|2"T)) . with the other factors accounting for 
the square root singularity of the BCS density of states. 
We now want to show that the identification 7 = is 
correct at next to leading order. To do so, we initially 
assume that the left /right gaps are different, so that the 
logarithmic divergence is absent and the perturbative ex- 
pansion in t is justified. Next, we keep only those terms 
that would become logarithmically divergent in the limit 
of equal gaps. 

To begin our derivation, we note that in Eq. (|A9[) the 
first term in square brackets multiplying Af lm originates 
from ((er+dj^d^)), as explained in Appendix El To- 
gether with the other term in square brackets, they give 
rise to the pure dephasing rate term in the master equa- 
tion (|A9[) via the equality 



2t 



n.rn.a 



\A d 



A d * 



= iT^ P+ {t) (B3) 



In what follow we first consider in some detail the next 
order contributions to ((er+d^d^)) and then discuss 
briefly the contributions to other averages. Without in- 
voking the lowest order Born approximation, the equa- 
tion of motion for ((<r + dj^d^)) is obtained by adding 
to the right hand side of Eq. (|A5[) the terms 



Appendix B: Dephasing at next-to-leading order 



The self-consistent equation ([28]) for requires go- 
ing beyond the lowest order (in the tunneling amplitude 



1 



A d -N a>p 



A d *M a ' p ■■ 



■ nm,ij 



Z A r * p a 'P 
2 W nm,ij 



2 ^ij ^nm,ij 



4- _ 4 r * R a i p 

9 i>3 nm,ij 



(B4) 
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with the definitions 

MZ,ij = P + {«f}&i P ^A})) (B5) 

- S ni S mj 6 apP+ [tf (1 - /«) + (1 - 

K£« = «<^ + *})) (B6) 

- «** {&j?fif P ; (B7) 

- S„iS mj S op p a [/„ L (1 - /*) + (1 - 

= «** {«f>fp; «ni«ma})> (B8) 



At lowest order, all the averages in the right hand side 
of Eq. (jBlljl vanish; non- vanishing contributions can in 
principle be found by considering once again the equa- 
tion of motions for those averages. As it is well known, 
proceeding in this manner we would obtain a hierarchy of 
coupled equations. — Here we make two approximations: 
first, we truncate the hierarchy at this level; second, as 
explained above we keep only those terms that in the 
limit of equal gaps would give logarithmically divergent 
contributions to the master equation. As a first step, this 
amounts to performing a mean-field like approximation 
in which the averages in the right hand side of Eq. (|Blip 
are written in terms of product of averages as in the fol- 
lowing example: 



R a ' p ■ ■ 

nm,ij 



a m a L -a L ^a R 

jp ipi Tier L mcr 



Sni^mjSap [f m f n ] 



q<r,P — 



a L] a R \a L ^a R 

ip 3 pi ri<7 ma 



(B9) 



(BIO) 



In introducing these definitions we have subtracted out 
the lowest order contributions already appearing in 
Eq. (|A5[) . Then in that equation and in Eqs. (|B5|) and 
(|B7[) the density matrix should be understood as the 
lowest (zeroth) order one. In other words, by construc- 
tion the quantities defined in Eqs. (|B5[) - (|B10[) account for 
higher order (in t) contributions; these can be found by 
considering the equations of motions for those quantities, 
such as 



id t M nm,ij = (^10 + 4. 



+ e?-ef)M£ 



k,l,p 



(Bll) 



l r 

2 



A kl a lfj, a kn> • M "nm,ij 



where A^f^ y stands for the anticommutator 



A a ' p 



{ & J 



m a L -a L ^a R 



(B12) 



L t/:,-R. A<r,P 



■A 



2((a+a^a R J((A P i%) 



where 



(B13) 



((A°' p 



7im,t3 l 



^ni^mj^ap \_fn 



/*) + (1 - fn )fn 



(B14) 

Similar expressions can be written for the other averages 
appearing in Eq. (jBllj) . In the second step we check 
which of the terms obtained in this way are logarith- 
mically divergent in the limit of equal gaps and discard 
those that are finite (here we employ again the Born- 
Markov^, and rotating wave approximations). 



Applying the above procedure to Eq. (|B11|) we find 
that the terms in the last two lines can be neglected, while 
in terms originating from the second line we use Eq. (|B3j) 
as well as Eq. (|A7[) (in the rotating wave approximation, 
we only need to keep the term in the right hand side of 
that equation that contains p+). Solving the equation 
for i ■ so obtained we finally arrive at 



= -«>+(<) SniSmjSrp [/„ L (1 - /*) + (1 - /„ L )/£] - 2t 2 p+(t)A%A^ m - f R ) + (1 - 

X [#(1 - frn) + (1 - fn)f*} e^~^f -ef+«>+)(t- U ) C dr (>£-^+«)+)(«-t) + ^ +ef -W0+){„-r)) 

(B15) 



We then use the same approach to find the expression 
for y [Eq. (|B6[) ]. which has the structure similar to 
that of the last term in Eq. (|B15j) . Using these results 



we get 



E( A d /V CT ' P 
\ ij nm,ij 

i,j,P 



A d * M a,p N 

ij nm.ij i 



= -I>+(t)^* m [tf (1 - f R ) + (1 - tf)f%\ (B16) 



t+ I dr e 

o 



i(e£-e£+»0+)(i— r) 
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To obtain the last term in curly brackets we used the 
identity 



/•t />U pt pt 

du dr h(r, u) — j dr j du h(r, u) 
Jo Jo 



JO JO 

pt PU 

— j du j dr h(u, r) 
Jo Jo 



(B17) 



to combine contributions coming from tj an d 

ij m a compact form. 
Using the same procedure one can find the expressions 
for the quantities defined in Eqs. (|B7 j) - (|B10j) . Those 
quantities, however, do not contribute to the master 
equation within the approximations we are employing (in 
particular, we remind that in the rotating wave approx- 
imation we neglect by assumption terms small by the 
factor r^/wio). Therefore, we obtain the following next- 
to-leading order equation of motion for ((a+ajj^a^.)) by 
substituting Eq. (|B16|) into Eq. (1B4[) and adding the re- 
sult to the left hand side of Eq. (|A5I) : 



- id t ({a+ = (wio + e L n ~ efn) 

- ~>C [(1 + P*) fn (1 - O (1 - P.)(l - /n L )/£] 

(B18) 

As explained at the beginning of this Appendix, we want 
to show that this equation agrees at next to leading or- 
der with the smearing obtained by introducing a finite 
decay rate in the terms responsible for dephasing, with 
the decay rate given by itself. Indeed, introducing 
this decay in Eq. (|A7[) we find 

Jo 

p+(r) - /*) + (1 - ft)f*\ 



[(1 + p z (r)) /„ L (1 - /*) - (1 - Pz (r))(l fk)f^ 



(B19) 

Taking the time derivative of this equation we get 

- id t {{<t + & n L a a* a )) = (u 10 + e L n - e£) ((a+ a^a^J 

+ iAt nP+ wo. - o + (i - .am 

- [(1 + P>) fn (1 - /™) - (! - P»)(l - /n )/£] 



(B20) 



T^A d n * m / drp + (r)e 4( " 10+e "- e ™ +ir * )( ^ T) 
Jo 

x [/„ L (i - /*) + (i - /M 



At next-to-leading order, one should expand the expo- 
nentially decaying part of p + [cf. Eqs. (f2"B")) -(f2"4" |) ] in the 
second line of Eq. (|B20|) and hence substitute there, with 
logarithmic accuracy, p + — > p+(l — T<pt). The last term 
in Eq. (|B20[) is explicitly of higher order, so one can use 
P+{t) ~ e 4Wl ° r and drop in the exponent. In this way 
we recover Eq. (IB18|) . thus showing for ((<7 + d^d^j CT )) the 
validity of our Ansatz. To complete the proof, we repeat 
the above steps for other averages, such as al^a^a)) 
and ((o^aJ^A^J). The latter contributes to the l/2Ti 
term in the master equation (|21[) and at next-to-leading 
order the only correction we find is that corresponding to 
the expansion of the exponentially decaying part of p+, 
as discussed above for the second line in Eq. (|B20I) . 



Appendix C: Phase relaxation in flux qubit and 
fluxonium 



1. Flux qubit 

In a flux qubit, the external flux threading the su- 
perconducting loop is tuned to half the flux quantum, 
/ = & e /&o — 1/2, and the potential energy takes the 
form of a double well. Then the qubit states |±) are 
the two lowest tunnel-split states in this potential with 
energy difference 



wio(/) = V e " 2 +[( 27r ) 2 ^(/- 1 /2)]' 



(CI) 



where for Ej 3> Ec we have 



e = 2^V^W } (S^j V4 e V^77^ (C2) 

Expressions for the renormalized parameters Ec, Ej, 
and El in terms of the bare parameters of the Hamilto- 
nian ([3]) can be found in Sec. IV. B of Ref.@. It was shown 
there that the matrix element A r s vanishes at / = 1/2 
because of symmetry considerations, thus leading to a 
minimum for the relaxation rate. Here we focus on the 
case / = 1/2 and therefore we need to evaluate the con- 
tribution to relaxation originating from the last line in 
Eq. (|16[) . The relevant matrix element is 



\Al\ = 



(C3) 



which equals unity at / = 1/2. Then from Eq. (fT5|) we 
obtain 




-E, 



neT P 



-A/T e 



(C4) 



Turning now to the dephasing rate, we find at / = 1/2 
the following expression for the matrix element 



\A. \ = 



D 



2V2 Ej \E C 



Ej 



1/3 



(C5) 
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where Da 1.45 is a numerical coefficient. 6 Using Eq. ([35 
and (|CJ4|) , after straightforward algebra we arrive at 



(C6) 



Due to the exponential suppression of the splitting, 
Eq. (IC2|) . this quantity is in general small. Indeed, 
for E c /A, T e /A > 0.01 and Ej/E c > 15 wc find 
2Tir,<, < 0.01. 




flux qubit one^ we find again that the pure dephasing 
rate is small compared to the relaxation rate for large 
Ej/Ec- The latter condition is not satisfied experimen- 
tally, since typically^ Ej/Ec < 5, and numerical cal- 
culations beyond the scope of the present work may be 
needed to address this parameter regime. However, we 
note that in all cases studied here decreasing the ratio 
Ej/Ec increases the relative contribution of pure de- 
phasing to 1 jT-i . 



Appendix D: Andreev bound states and ionization 

rate 



2. Fluxonium 

In the fluxonium an array of Af > 1 identical junc- 
tions (each with Josephson energy Eji ^> Eci large 
compared to their charging energy) acts as a lossy in- 
ductor connected to a weaker junction with Ejq < Ej\. 
The inductive energy of the array is El — Eji/M and 
the losses are due to quasiparticle tunneling through the 
array junctions. In fact, for external flux near half the 
flux quantum the relaxation time is determined by this 
loss mechanism^ 

a/t b ( mo (1/2) y (c) 
wio(/r v wio(/) J ' [ ' 

since as discussed above for the flux qubit the contribu- 
tion of the weaker junction is suppressed at / = 1/2 
[cf. Eq. |(C4|)]. Note that at / = 1/2 the rate in 
Eq. (IC7l) is larger than that in Eq. (|C4)) by the factor 
(A/Ej)(E l /cj 10 (1/2)). 

To calculate the dephasing rate, we note that at / = 
1/2 the matrix element for the weak junction is the same 
as for the flux qubit, Eq. (|C5I) . 




The goals of this Appendix are to derive Eqs. (|5"T)l 
starting from the model defined by Eqs. (H)-©, an d to 
estimate the ionization rates due to qubit-quasiparticlcs 
interaction and flux noise. In the low energy limit where 
the characteristic energy of the quasi-particles SE as well 
as the qubit transition frequency uio are small compared 
to the superconducting gap A, we approximate the BCS 
coherence factors as u J n w v 3 n m l/y/2. Then considering 
for now a single channel junction, Eq. © takes the form 6 



Ht = it sin(<£/2) 



a L ^a R 

ncr mcr 



H.c. 



(Dl) 



Assuming for simplicity identical left/right leads, we per- 
form a canonical rotation into a new quasiparticle basis 
defined by the operators 



7±r. 



1 

V2 



(D2) 



In this basis we have [cf. Eq. 



H qp — H qp+ + H qp _ , H qp ± - y2 e n 7±„ CT 7±no- (D3) 



n.(T 



Ht = t sin 



-no-7— ma J+nuJ+mo- 



(D4) 



, d , = D mo(l/2) (Ejq 
1 sM 2V2 Ejo \E C0 

while for each array junction we get 



1/3 



I A d I 



2M 



(C8) 



(C9) 



Then the coefficient containing the sum over all junctions 
is 



M 
3=0 



Ejq + ^El) \Ai 



«,0 ' 



(C10) 



which in the limit Ej/Ec 3> 1 is exponentially sup- 
pressed, see Eq. (|C2|) . Therefore at / = 1/2 the de- 
phasing rate has the same exponential suppression 
as in the flux qubit. Since as discussed above the flux- 
onium relaxation rate is parametrically larger than the 



Denoting with \j) and £j the eigenstates and eigenener- 
gies of H(p [Eq. (jjJJ)], the total Hamiltonian H can then be 
split into parts that are respectively diagonal and non- 
diagonal in the qubit subspace, H = + H n d, with the 
diagonal part defined as 

H d = £ Sj \j) U> U I (D5) 



where 

H 3 ± = H qp± Ttsjj 7± na J±ma- (D6) 

and we have used the definition (I10p for the matrix ele- 
ments Sij . The non-diagonal part is given by 



-rial— ma 1+nal+ma 

(D7) 
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It describes real transitions in which qubit and quasipar- 
ticles exchange energy. The term proportional to t in the 
diagonal part, on the other hand, accounts for virtual 
transitions that renormalize the spectrum. Indeed, as we 
show next, for Sjj > (sjj < 0) there exists a sub-gap 
Andreev bound state in the 7+ (7-) subspace. Because 
the two subspaces are uncoupled, we can restrict our- 
selves to either one of those; in the following we consider 
the 7+ subspace. 

Since is diagonal in the qubit space, to find the 
spectrum we only need to calculate the eigenvalues of 
the quasiparticle Hamiltonians Hj±. We denote with 
\Aj) the wavefunction of the Andreev state when the 
qubit is in state \j); to solve the Schrddinger equa- 



tion Hj+\Aj 



E\Aj) we make the Ansatz \Aj 



;ncr "'jn i+na\^)-> where |0) denotes the quasiparticle vac- 
uum state, 7± nC r|0) = 0, and obtain the following system 
of linear equations 



ts 



1 



33 



E 



(D8) 



To find the eigenenergy E, we sum both sides over n and 
in the low energy limit we write e n s» A + f„ 2 /(2A); 
then the sum over n in the right hand side can be ap- 
proximated by an integral, ^ ra f=a vq J (if, and we arrive 
at 



1 = TtVQts 



33 



2A 



E 



(D9) 



A solution with energy E < A exists if and only if Sjj > 
(the opposite holds in the 7_ subspace where a bound 
state exists if and only if Sjj < 0.). The corresponding 
bound state energy is 



Ef = A[l — 2(TTVot) 2 Sjj 2 ], 



(D10) 



This energy depends on the state of the qubit via the 
matrix element Sjj. However, for the low-energy states 
of the phase qubit this matrix element is the same at 
leading order in Ec/^>w *C 1, since the square of the 
matrix clement is£ 



1 



} Ec 
Ec r x 

— bkj-i 



2 



+ (j + l)Si,j + l] cos 2 ^ 



(Dll) 



up to higher order terms oc (Ec /wio) 2 . Keeping only the 
leading term in this equation, introducing the transmis- 
sion probability T = (27Wot) 2 in Eq. (ID10I) . and general- 
izing it to multiple channels, we arrive at Eq. (|57| . (In 
that equation the subscript p denotes the transmission 
channel, and we have dropped the qubit state index j 
since, as explained above, the leading order expression is 
independent of j.) 

For later use, we note that the normalization condi- 
tion J2n( a jn) 2 ~ 1/2, which accounts for spin degener- 
acy, together with the square of Eq. (|D8I) . leads to the 



amplitudes 



1 (2At^) 3 / 4 



3 n 



2Aw 



A ' 



(D12) 



where ujf — A — E A is the binding energy. 



1. Ionization rate 

The ionization of the Andreev level can be caused by 
quantum fluctuations of the phase difference across the 
junction induced by the finite charging energy Ec", the 
ionization rate can be calculated using Fermi's golden 
rule by treating the non-diagonal part (ID7[) of the Hamil- 
tonian as a perturbation. For a qubit initially in the state 
the ionization rate T A is given by 



Tf = 2 7 r^|<i,e in |ff nd |i,A i 



Here 



x S(£ 3 Ef)(l - f(e n )). (D13) 

is a scattering state in the continuum part 
of the quasiparticle spectrum and the factor (1 — 
gives the probability that this state is empty. The matrix 
element in Eq. (|D13[) is the product of the off-diagonal 
matrix element sji times the overlap of the wavefunctions 
of bound and scattering states at the junction, 



(j,e jn \H nd \i,Ai) = 



-isu «/>!L (e n ) as. 



(D14) 



where 4>jm(^n) — (^m\^jn) and \e m ) are the eigenstates 
of H qp+ , see Eq. (|E)3|) . Next we calculate the wavefunc- 
tions for the continuum states by solving the scattering 
problem in the standard T-matrix approach^ 

We focus again on the 7+ subspace and write Hj+ = 



<qp+ +-H3I With Hj! = -tSjj Enm.a 7+na7+m<7 [cf. 



Eq. (|D6|) ] . From the Schrddinger equation, we have for 
the scattering states \cj n ) 



1 



■Hij\ejn 



f i0 
iQ + 



\e n ), (D15) 



where we have defined the T-matrix as 



T? ( e «) = ^1? + Hij [e n — H< 



(Die) 



The T-matrix is related to the quasiparticle Green's func- 
tion G, via 



Gj = 9 + 9^j9 , 



(D17) 



where g is the (diagonal in momentum) bare quasiparticle 
Green's function g n (u>) — — e n + i0 + ). Using the 
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inverse of this equation: Tj = g 1 Gjg 1 — g , we find 
upon projecting Eq. (|D15|) onto \e m ) 



fpjmi^n) = lim Gj. mn (uj)(g n {uj)) 1 . 



(D18) 



The Green's function, as obtained from the equations of 
motion for 7+n CT , is given by 



Gj,nm{k>) — 5 nm g n (ui) — — 



isjjg n (oj)g m (oj) 



Hence the continuum states are 



1 — mVQtSjj 



2A 



(D19) 



(D20) 



where we have used that in the low energy limit 
J2 P 9p( e 7i) ~ -iirv a v / 2A /(e n - A). 
Using Eqs. (|D12|) and (|D20|) we find 



Defining the frequency shift 8ui q = Uxo — oj q and using 



Eq. ()59|) to estimate its value, we may eliminate N a , 
and rewrite the above as 



Tir 



A 
tot 



a A/T e 



Ej Suj 



N c hT e uj w 



^1 „ 4 x i 4^i ( D26 ) 



where we used Ej/N ch ~ 1(T 5 A and T e 140 mK. Thus 
when StUq > 10~ 4 wio, the qubit relaxation is likely dom- 
inated by the ionization process, rather than by quasi- 
particle transitions within the continuum. However, we 
note that the typical shift is much smaller than this, 



Sujq/uJia ~ e 



-A/T e 



3 x 1(T 7 , i.e. p w 0.012, so the 



contribution of ionization to qubit relaxation is negligi- 
ble unless N occ is anomalously large. 



2. Ionization by flux noise 



y^Qw = 2 1/4 7w \/iAs ii , 



1 + iTrvotSjj x /2A/(e n - A) 



(D21) 



(D22) 



Substituting these expressions into Eq. (|D13|) . and con- 
sidering explicitly the case of a phase qubit, using the ex- 
pressions for the matrix elements given in Eq. (|D11[) 
finally yields for the ionization rate of a single-channel 
junction 



= J^-M 1 +cos tfo)- 
8oj w 



2ur. 



xil-m-Sj-x+Ef)), (D23) 



with tpo, uip, and luiq defined in Sec. IIV Al and we used 
that for a single-channel junction Ej — A^vot) 2 . 

The above result can be easily generalized to the case 
of N c h independent channels. Assuming for simplicity 
identical transmission amplitudes, the Andreev binding 
energy can be written as co A = 2EjSjj 2 /N c h which for a 
phase qubit reduces approximately to u> A « Ej/N c h- We 
assume N c h ^> 1 sufficiently large so that uj a <C u>w and 
obtain for the ionization rate of each occupied channel 



1 



u! p 2 / 2uj a 1 + cos ipo 



4:N ch wio V w io 



(D24) 



A single ionization event is sufficient to relax the 
qubit energy, and the probability of at least one ion- 
ization event taking place during time t, when initially 
N occ < N c h Andreev levels are occupied, is given by 
p = 1 — e~ Nocari *. Introducing the total ionization rate 
A and using Eq. (1551) . leads to the estimate 



L tot 



NoccTi 



TiT A t 




(D25) 



As an example of an extrinsic ionization mechanism, 
we consider here low frequency (<C wio) flux noise. Small 
fluctuations <5<& e (t) <C $o of the external flux induce 
small fluctuations <P\{t) of the phase difference ipo, 



pi(t) = 2tt 



E L 



$0 E L + Ej COS tpo 



(D27) 



[see Eq. (|36p ]. Since the low- frequency fluctuations do 
not induce qubit transitions, their effect is accounted for 
by substituting ipg — > ipo + <p>i(t) into the diagonal ma- 
trix element Sjj in Eq. (|D6|) . At linear order in ipi we 
thus obtain the time-dependent perturbation (in the 7+ 
subspace) 



V-(i) = -<^co 8 



(f) 



^ y 7+ricr 



7+r, 



(D28) 



Using Fermi's golden rule and following similar steps as 
in the previous section, the total ionization rate can be 
expressed as 



1 tot ~ iv occ 



El 



ch 



3/2 



sin — 
2 



1 + COS ifo 



(D29) 



x / duo S w {u)^-^{1 - f{u + E A )). 



UJ 



where S vv (oj) = 1/(2tt) / e^^i^^O)}^ is the phase 
fluctuation spectrum and the binding energy introduces a 
natural low-frequency cutoff. For non-degenerate quasi- 
particles, f(u>A + Ea) <C 1, and a power-law spectrum of 
= (Sip) 2 1 \2iruj a ) , we obtain 



2tt 



3/2 



1 



COS tpo 



dx 



^/x~\ 



(D30) 



17 



For a = 1 (pure 1// noise) 
equal to 7r/2 and since uja = 
arrive at 



the remaining integral is 
2(Ej/N ch )sin 2 (<f /2), we 



rA (M 2 



N . 



Ej \ 1 + cos (f 



ch 



(D31) 



The measured 2 ^ magnitude of the fluctuations is small, 
~ 10" 6 ; since Ej/N ch ~ 1CT 5 A and N occ < AT c/l < 
10 7 , we estimate this rate to be much smaller than 1 Hz. 



Appendix E: Modifications of the density of states 

The logarithmic divergence of the daphasing rate and 
its regularization discussed in Sec. Mil are a consequence 
of the square root singularity of the BCS density of states 
at the gap edge. Here we discuss two other mechanisms 
that also can regularize the divergence and show that for 
Al-base qubits used at present they do not modify the 
results in the main text. 

To begin with we consider the broadened density of 
states introduced by Dynesi^ to interpret experimental 
tunneling data. This phenomenological density of states 
is characterized by a broadening parameter <C A and 
a finite density of subgap states. These states give rise to 
an additional contribution to the dephasing rate which we 
denote with assuming quasi-equilibrium, it is given 

20 



by^i 



I7(T e ) 



A 



Te 

A 



(El) 



and it is always smaller than the broadening, T s ^ <C To- 
Comparing Eqs. and (|E1[) . we see that a small broad- 
ening in the latter can compensate for the exponential 
suppression of the quasiparticle occupation in the former. 
Then we can distinguish three regimes: 1. at "high" 
temperatures, the dephasing rate is given by Eq. ([35]) . 
since the broadening can be neglected in calculating T^. 
The high-temperature regime is defined by the condition 
Td ^ ^4>(T e ); 2. at intermediate temperatures, when 
T T( T e) < r (T e ) < T D , the broadening of the density 
of states cannot be neglected. With logarithmic accu- 
racy, this amount to substitute — > To in the last 
term in Eq. (|33[) [and hence replace the square bracket 
in Eq. with ln(T e /r£>); we note that since this sub- 
stitution affects only the logarithm, use of Eq. (|35|) still 



gives a correct order-of-magnitude estimate]. 3. at low 
temperatures, such that Y,p(T e ) < T s ^(T e ) the subgap 
contribution becomes dominant. 

In recent measurements^ the intrinsic value of the 
broadening parameter in aluminum was found to be 
small, Try/ A < 2 x 1CP 7 . Using this value and the re- 
sults of the next section, our estimates show that the 
low-temperature regime is entered for T e < 60 mK. In 
experiments with superconducting resonators^ 7 - as well 
as qubits 7 -^ the quasiparticle effective temperature is 
larger, T e ~ 140 mK, so we can neglect the subgap contri- 
bution to the dephasing rate for Al-based qubits, which 
we focus on in this paper. However, the subgap contri- 
bution may be relevant in other systems, such as qubits 
fabricated with niobium^ 

While the above consideration are based on a phe- 
nomenological model, an intrinsic modification of the 
continuum part of the density of states near the junc- 
tion is due to the presence of Andreev bound states. 
They modify the square root singularity into a square 
root threshold, 



2A 



/2AVw - A 

u - E A 



(E2) 



with E A the energy of the bound state defined in Eq. (|57|) 
(here we consider for simplicity the single channel case). 
The above substitution can be obtained using Eq. (|D19|) 
for the Green's function to calculate the density of states. 
Assuming the binding energy io A = A — E A to be small 
compared to the typical quasiparticle energy, oj a <C SE, 
we find that the substitution (|E2j) would lead to the 
replacement of with u> A in the right hand side of 
Eq. (|33|) . In quasi-equilibrium this amount to replacing 
the square brackets in Eq. (1331) with 



In- 



T, 



In. 



Ik 

Ej 



■hi Nr.. 



(E3) 



where 3> 1 is the number of channels in the junc- 
tion. We note that the tunneling limit we are considering 
consists in taking the transmission amplitude t — > at 
finite Ej, which implies N c h — > oo. Then in this limit 
the self-consistent approach is justified with logarithmic 
accuracy as explained in Appendix [B] . 
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